clear; clc;
% 在不同 AE (0.02, 4.4) 下，能量 H(t) 的时间序列，平均能量 <H>.

% 数设定
a     = 0.7; 
b     = 0.8; 
c     = 0.1; 
alpha = 0.25; 
beta  = 0.01; 
delta = 0.1;

% 初始条件
y0 = [0.2; 0; 0.01];

% 积分时间
T_total = 1200;   % 总模拟时间
T_trans = 400;    % 舍去瞬态时间，后 800 时间单位用于计算 <H>
h       = 0.001;  % RK4 步长(可根据需要调整)

% AE 分别是 0.02 和 4.4
AE_list = [0.02, 4.4];

tlo = tiledlayout(2, 1);

for i = 1:2
    AE_val = AE_list(i);

    % 将参数打包
    params = struct('a', a, 'b', b, 'c', c, 'alpha', alpha, ...
                    'beta', beta, 'delta', delta, 'AE', AE_val);

    % 使用 hybridNeuron_AE 数值积分
    [t, Y] = rk4(@(tt,yy) hybridNeuron_AE(tt, yy, params), [0, T_total], y0, h);

    % 提取 x, y, z
    x_data = Y(:,1);
    y_data = Y(:,2);
    z_data = Y(:,3);

    % 论文公式 (7) 计算能量 H(t) - 
    % H = 0.5*x^2 + (1/(2c))*y^2 + 0.5*b*z*sin(beta*z)*y
    H_data = 0.5 * x_data.^2 ...
           + (1/(2*c)) * (y_data.^2) ...
           + 0.5 * b .* z_data .* sin(beta .* z_data) .* y_data;

    % 舍去瞬态 (t > T_trans)，计算后期平均能量
    idx = (t > T_trans);
    t_plot = t(idx);
    H_plot = H_data(idx);

    mean_H = mean(H_plot);

    nexttile(i);
    plot(t_plot, H_plot, 'Color', [0, 0.6, 0], 'LineWidth', 1.5);  % 绿色曲线
    xlabel('\tau');
    ylabel('H');
    if i == 1
        title(sprintf('(a) A_E=%.2f, <H>=%.5f', AE_val, mean_H));
    else
        title(sprintf('(b) A_E=%.2f, <H>=%.5f', AE_val, mean_H));
    end
    grid on;
end
